## Code to replicate Table 1 of Abrajano, Elmendorf, and Quinn. "Measuring
## Perceived Skin Color: Spillover Effects and Likert-type Scales". JOP.


library(clubSandwich)


sink(file="AbrajanoElmendorfQuinn-Table1-output.txt")

## Part 1
## E[y | d = W] > E[y | d = B]
##
## where the expectation is over all respondents assigned to M1 or M2,
## all photos, and all target positions

mydata <- read.csv("ScaleRaceSpring2017clean.csv")

## just black distractor photos
mydata.b <- mydata[mydata$condition %in% c("P1B", "P2B"),]
## just white distractor photos
mydata.w <- mydata[mydata$condition %in% c("P1W", "P2W"),]

## keep must the target photos
target.inds <- c(6, 7, 10, 12, 15, 17, 19, 20, 22, 23)
keep.vars <- c(paste("X", target.inds, "_Q339", sep=""),
               paste("X", target.inds, "_Q329", sep=""))

mydata.b <- mydata.b[, keep.vars]
mydata.w <- mydata.w[, keep.vars]

## reshape into long format
mydata.b.long <- reshape(mydata.b, direction="long", varying=1:20, v.names="MM")
mydata.w.long <- reshape(mydata.w, direction="long", varying=1:20, v.names="MM")

## keep only obs with non-missing values for MM
mydata.b.long <- na.omit(mydata.b.long)
mydata.w.long <- na.omit(mydata.w.long)

lm.b.out <- lm(MM ~ 1, data=mydata.b.long)
lm.w.out <- lm(MM ~ 1, data=mydata.w.long)

tab.b <- coef_test(lm.b.out, vcov="CR2", cluster=mydata.b.long$id,
                   test="Satterthwaite")

tab.w <- coef_test(lm.w.out, vcov="CR2", cluster=mydata.w.long$id,
                   test="Satterthwaite")
WminusB.diff <- tab.w[1,1] - tab.b[1,1]
WminusB.var <- tab.w[1,2]^2 + tab.b[1,2]^2
WminusB.SE <- sqrt(WminusB.var)
WminusB.z <- WminusB.diff / WminusB.SE
WminusB.pval <- 2*(1-pnorm(abs(WminusB.z)))

cat("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
cat("Part 1:  E[y | d = W] - E[y | d = B] > 0\n")
cat("\n")
cat("----------------------------------------------------------------\n")
cat("Avg. MM score among respondents with white distractors\n")
print(as.data.frame(tab.w), digits=4)
cat("\nAvg. MM score among respondents with black distractors\n")
print(as.data.frame(tab.b), digits=4)
cat("\n\n----------------------------------------------------------------\n")
cat("WminusB.diff     WminusB.SE     WminusB.z    p-value (2-tailed)\n")      
cat("----------------------------------------------------------------\n")
cat(WminusB.diff, "       ", WminusB.SE, "    ", WminusB.z, "       ",
    WminusB.pval, "\n")
cat("----------------------------------------------------------------\n")
cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n\n")











## Part 2
## E[y_s | d = W] > E[y_s | d = B]
## where the expectation is over all respondents assigned to M1 or M2,
## all photos, and conditional on s = 6.

mydata <- read.csv("ScaleRaceSpring2017clean.csv")

## just black distractor photos
mydata.b <- mydata[mydata$condition %in% c("P1B", "P2B"),]
## just white distractor photos
mydata.w <- mydata[mydata$condition %in% c("P1W", "P2W"),]

## keep just the target photos
target.inds <- c(6)
keep.vars <- c(paste("X", target.inds, "_Q339", sep=""),
               paste("X", target.inds, "_Q329", sep=""))

mydata.b <- mydata.b[, keep.vars]
mydata.w <- mydata.w[, keep.vars]

## reshape into long format
mydata.b.long <- reshape(mydata.b, direction="long", varying=1:2, v.names="MM")
mydata.w.long <- reshape(mydata.w, direction="long", varying=1:2, v.names="MM")

## keep only obs with non-missing values for MM
mydata.b.long <- na.omit(mydata.b.long)
mydata.w.long <- na.omit(mydata.w.long)

lm.b.out <- lm(MM ~ 1, data=mydata.b.long)
lm.w.out <- lm(MM ~ 1, data=mydata.w.long)

tab.b <- coef_test(lm.b.out, vcov="CR2", cluster=mydata.b.long$id,
                   test="Satterthwaite")

tab.w <- coef_test(lm.w.out, vcov="CR2", cluster=mydata.w.long$id,
                   test="Satterthwaite")
WminusB.diff <- tab.w[1,1] - tab.b[1,1]
WminusB.var <- tab.w[1,2]^2 + tab.b[1,2]^2
WminusB.SE <- sqrt(WminusB.var)
WminusB.z <- WminusB.diff / WminusB.SE
WminusB.pval <- 2*(1-pnorm(abs(WminusB.z)))

cat("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
cat("Part 2:  E[y_6 | d = W] - E[y_6 | d = B] > 0\n")
cat("\n")
cat("----------------------------------------------------------------\n")
cat("Avg. MM score among respondents with white distractors\n")
print(as.data.frame(tab.w), digits=4)
cat("\nAvg. MM score among respondents with black distractors\n")
print(as.data.frame(tab.b), digits=4)
cat("\n\n----------------------------------------------------------------\n")
cat("WminusB.diff     WminusB.SE     WminusB.z    p-value (2-tailed)\n")      
cat("----------------------------------------------------------------\n")
cat(WminusB.diff, "       ", WminusB.SE, "    ", WminusB.z, "       ",
    WminusB.pval, "\n")
cat("----------------------------------------------------------------\n")
cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n\n")








## Part 3
## E [y | d = W, i in I_M1 ]  -  E [y | d = B, i in I_M1 ] >
##  E[y | d = W, i in I_M2 ]  -  E [y | d = B, i in I_M2 ] 
## where the expectations are over all respondents assigned to M1 or M2,
## all photos, and all target positions.

mydata <- read.csv("ScaleRaceSpring2017clean.csv")

## just black distractor photos
mydata.b <- mydata[mydata$condition %in% c("P1B", "P2B"),]
## just white distractor photos
mydata.w <- mydata[mydata$condition %in% c("P1W", "P2W"),]

## keep just the target photos
target.inds <- c(6, 7, 10, 12, 15, 17, 19, 20, 22, 23)
keep.vars.template <- paste("X", target.inds, "_Q329", sep="")
keep.vars.notemplate <- paste("X", target.inds, "_Q339", sep="")

mydata.b.m1 <- mydata.b[, keep.vars.template]
mydata.b.m2 <- mydata.b[, keep.vars.notemplate]
mydata.w.m1 <- mydata.w[, keep.vars.template]
mydata.w.m2 <- mydata.w[, keep.vars.notemplate]

## reshape into long format
mydata.b.m1.long <- reshape(mydata.b.m1, direction="long",
                            varying=1:10, v.names="MM")
mydata.b.m2.long <- reshape(mydata.b.m2, direction="long",
                            varying=1:10, v.names="MM")
mydata.w.m1.long <- reshape(mydata.w.m1, direction="long",
                            varying=1:10, v.names="MM")
mydata.w.m2.long <- reshape(mydata.w.m2, direction="long",
                            varying=1:10, v.names="MM")

## keep only obs with non-missing values for MM
mydata.b.m1.long <- na.omit(mydata.b.m1.long)
mydata.b.m2.long <- na.omit(mydata.b.m2.long)
mydata.w.m1.long <- na.omit(mydata.w.m1.long)
mydata.w.m2.long <- na.omit(mydata.w.m2.long)

lm.b.m1.out <- lm(MM ~ 1, data=mydata.b.m1.long)
lm.b.m2.out <- lm(MM ~ 1, data=mydata.b.m2.long)
lm.w.m1.out <- lm(MM ~ 1, data=mydata.w.m1.long)
lm.w.m2.out <- lm(MM ~ 1, data=mydata.w.m2.long)

tab.b.m1 <- coef_test(lm.b.m1.out, vcov="CR2", cluster=mydata.b.m1.long$id,
                   test="Satterthwaite")
tab.b.m2 <- coef_test(lm.b.m2.out, vcov="CR2", cluster=mydata.b.m2.long$id,
                   test="Satterthwaite")

tab.w.m1 <- coef_test(lm.w.m1.out, vcov="CR2", cluster=mydata.w.m1.long$id,
                   test="Satterthwaite")
tab.w.m2 <- coef_test(lm.w.m2.out, vcov="CR2", cluster=mydata.w.m2.long$id,
                   test="Satterthwaite")

WminusB.m1.diff <- tab.w.m1[1,1] - tab.b.m1[1,1]
WminusB.m1.var <- tab.w.m1[1,2]^2 + tab.b.m1[1,2]^2
WminusB.m1.SE <- sqrt(WminusB.m1.var)
WminusB.m1.z <- WminusB.m1.diff / WminusB.m1.SE
WminusB.m1.pval <- 2*(1-pnorm(abs(WminusB.m1.z)))

WminusB.m2.diff <- tab.w.m2[1,1] - tab.b.m2[1,1]
WminusB.m2.var <- tab.w.m2[1,2]^2 + tab.b.m2[1,2]^2
WminusB.m2.SE <- sqrt(WminusB.m2.var)
WminusB.m2.z <- WminusB.m2.diff / WminusB.m2.SE
WminusB.m2.pval <- 2*(1-pnorm(abs(WminusB.m2.z)))

WminusB.diff.diff <- WminusB.m1.diff - WminusB.m2.diff
WminusB.diff.diff.var <- WminusB.m1.var + WminusB.m2.var
WminusB.diff.diff.SE <- sqrt(WminusB.diff.diff.var)
WminusB.diff.diff.z <- WminusB.diff.diff / WminusB.diff.diff.SE
WminusB.diff.diff.pval <- 2*(1-pnorm(abs(WminusB.diff.diff.z)))


cat("\n@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n")
cat("Part 3: (E [y | d = W, i in I_M1 ]  -  E [y | d = B, i in I_M1 ]) -\n")
cat("         (E[y | d = W, i in I_M2 ]  -  E [y | d = B, i in I_M2 ]) > 0\n") 
cat("\n")
cat("----------------------------------------------------------------\n")
cat("Avg. MM score among respondents with white distractors & no template\n")
print(as.data.frame(tab.w.m1), digits=4)
cat("\nAvg. MM score among respondents with black distractors & no template\n")
print(as.data.frame(tab.b.m1), digits=4)
cat("\n\n")
cat("Avg. MM score among respondents with white distractors & template\n")
print(as.data.frame(tab.w.m2), digits=4)
cat("\nAvg. MM score among respondents with black distractors & template\n")
print(as.data.frame(tab.b.m2), digits=4)

cat("\n\n----------------------------------------------------------------\n")
cat("WminusB.m1.diff  WminusB.m1.SE  WminusB.m1.z   p-value (2-tailed)\n")      
cat("----------------------------------------------------------------\n")
cat(WminusB.m1.diff, "       ", WminusB.m1.SE, "    ", WminusB.m1.z, "       ",
    WminusB.m1.pval, "\n")
cat("----------------------------------------------------------------\n")

cat("\n\n----------------------------------------------------------------\n")
cat("WminusB.m2.diff  WminusB.m2.SE  WminusB.m2.z   p-value (2-tailed)\n")      
cat("----------------------------------------------------------------\n")
cat(WminusB.m2.diff, "       ", WminusB.m2.SE, "    ", WminusB.m2.z, "       ",
    WminusB.m2.pval, "\n")
cat("----------------------------------------------------------------\n")

cat("\n\n****************************************************************\n")

cat("----------------------------------------------------------------\n")
cat("WminusB.dif.dif  WminusB.dif.dif  WminusB.dif.dif.z   p-value (2-tailed)\n")      
cat("----------------------------------------------------------------\n")
cat(WminusB.diff.diff, "       ", WminusB.diff.diff.SE, "    ", WminusB.diff.diff.z, "       ",
    WminusB.diff.diff.pval, "\n")
cat("----------------------------------------------------------------\n")

cat("@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@\n\n\n")



sink()
